Dynamics analysis and optimal control of SIVR epidemic model with incomplete immunity

In this paper, we establish an SIVR model with diffusion, spatially heterogeneous, latent infection, and incomplete immunity in the Neumann boundary condition. Firstly, the threshold dynamic behavior of the model is proved by using the operator semigroup method, the well-posedness of the solution and the basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\Re _{0}$\end{document}ℜ0 are given. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\Re _{0}<1$\end{document}ℜ0<1, the disease-free equilibrium is globally asymptotically stable, the disease will be extinct; when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\Re _{0}>1$\end{document}ℜ0>1, the epidemic equilibrium is globally asymptotically stable, the disease will persist with probability one. Then, we introduce the patient’s treatment into the system as the control parameter, and the optimal control of the system is discussed by applying the Hamiltonian function and the adjoint equation. Finally, the theoretical results are verified by numerical simulation.


Introduction
The SARS in 2003, the Zika virus (ZIKV) invasion in 2013, the H7N9 avian influenza pandemic, and the emergence of the Dengue virus in the world, these recurrent infectious diseases and various emerging infectious diseases have been challenging modern life and medical standards [1]. For example, COVID-19, which broke out in 2019, is still affecting the world. As of August 24, 2021, the cumulative number of COVID-19 cases and deaths has reached 212,357,898 and 4,439,843. Therefore, how to prevent and control the occurrence and spread of infectious diseases is one of the hot issues today.
From the perspective of mathematics, the study of infectious diseases usually starts according to the transmission mechanism of diseases, which is analyzed by establishing mathematical models. The earliest epidemic model was established by Kermack and Mckendrick. They established the plague susceptibility infection removal model (SIR) [2] in 1927 and the plague susceptible infected susceptible model (SIS) [3] in 1932, respectively. Since the establishment of SIR and SIS models, most of the subsequent research is based on the standard SIR model. Among the existing prevention and treatment methods for infectious diseases, vaccine injection is one of the fast and effective methods. For example, in the prevention and control of COVID-19, vaccine injection can reduce the infection rate of the Delta variant virus to a certain extent. Therefore, an increasing number of researchers take vaccine injection into account in the process of modeling infectious diseases to make the model close to the actual situation.
Among them, Chen et al. in [4] established the susceptibility vaccination-infection isolation-recovery (SVIQR) model and susceptibility-vaccination-infection isolation (SVIQS) model, respectively. The basic reproduction number of the two models was given. Furthermore, the global attractivity and the global asymptotic stability of the solutions were proved by the Lyapunov function method. And the existence of backward bifurcation was also proved. In [5], Kribs-Zaleta and Velasco-Hernández studied a simple SIV model with inoculation and demonstrated the backward bifurcation of solutions to some parameter values. At the same time, complete bifurcation analysis of the model was given under the condition that the vaccine reduces the basic reproduction number. Liu et al. in [6] studied the following SIVR system: Here, λ is the constant update rate of the susceptible host and α is the rate at which susceptible individuals are vaccinated. β 1 and β 2 are the transmission rates of infected people in contact with susceptible groups and vaccinated groups, respectively. Since vaccinated individuals may have partial immunity during vaccination, it is assumed that β 1 > β 2 . μ is the host mortality per compartment. γ and δ are the recurrence rates of vaccinated people and infected people, respectively. All of these parameters are assumed to be positive. In [6], the authors gave the threshold dynamics of system (1.1) by using the basic reproduction number, showing that reducing the number of infected individuals by vaccination can control the disease. In addition, many researchers have studied infectious disease models with immunization from the perspective of age structure [7] and pulse vaccination [8].
The above models are all established in a homogeneous space environment. However, in practice, the transmission of some diseases is often related to spatial location. For example, the transmission rate of COVID-19 in Asia is different from that in North America. In [9], Wu et al. discussed a class of spatially heterogeneous host-pathogen models. The authors used the basic reproduction number to discuss the threshold dynamic behavior of the models and gave the asymptotic behavior of the models. In [10], a reaction-diffusion model of SVIR infection in a spatially heterogeneous environment was proposed. The authors gave the proof of the extinction and persistence of the disease by giving a basic reproduction number. In [11], the authors established an SIVS epidemic model with a degree-dependent transmission rate and incomplete vaccination on a scale-free network. The global asymptotic stability of the equilibrium and the global attractivity of the unique endemic equilibrium were proved. In addition, the effects of various immunization programs such as unified immunization, target immunization, and acquaintance immunization were studied and compared.
Motivated by the recent development of epidemic modeling, the optimal control problem is often discussed in some cases, optimal control theory is one of the important branches of mathematical optimization, which is often used to study how to find a control for a dynamic system in a period of time to optimize the objective function. Thus we consider two different models based on (1.1). The first is a direct extension of (1.1). A reaction-diffusion SIVR model is established based on spatial heterogeneity with incomplete immunity. The well-posedness of the system is discussed by using the operator semigroup method. At the same time, the global dynamic behavior of the system solution is discussed by analyzing the basic reproduction number. In the second model, as [12,13], it is assumed that the spread of a disease can be influenced by decision-makers. That is, decision-makers can control the response rate to a certain extent by increasing the treatment ability or the efficiency of drug treatment. Therefore, by further expanding the model, we obtain the control system under the assumption of limited control resources. Considering the progress of medical technology, the targeted treatment for patients will be gradually developed, so we will consider the targeted treatment for patients as a control parameter in the system and discuss its optimal control problem. In addition, we analyze the optimal control of the system by using the Hamilton equation and the adjoint equation.
In the proof, we may encounter the following problems: • How to determine the basic regeneration number of the system. In the presence of the diffusion term, it is necessary to select an appropriate method to represent the basic reproduction number 0 , which is an important prerequisite for discussing the dynamic behavior of the system by using 0 as the threshold value. • Can the existence of optimal control be obtained? Because of the existence of diffusion terms, it is difficult to define the adjoint equation and the Hamiltonian function of the control system. At the same time, there are some requirements for the selection of parameters in the numerical simulation. In view of the above problems, this article is organized as follows. In Sect. 2, an SIVR model with incomplete immunity and spatial heterogeneity is established. Furthermore, the well-posedness of the model is derived, meanwhile, the global existence and global attractiveness of the solution are proved. Section 3 is devoted to studying the threshold dynamic behavior of the system. The extinction or persistence of diseases is analyzed by using the basic reproductive number as the threshold. In Sect. 4, the optimal control of the control system is analyzed by taking the treatment for the patient as the control parameter in the system. Meanwhile, the optimal control problem is discussed by using the Hamiltonian function and adjoint equation. Finally, in Sect. 5, the corresponding results are verified by numerical simulation.

Model formulation and well-posedness
In this paper, the spatial heterogeneity of the spread of infectious diseases and spatial diffusion is considered. In addition, for vaccines, we consider vaccination rates in susceptible individuals and the effectiveness of the vaccine. Based on model (1.1), we can establish the following epidemic model of SVIR with incomplete immunity. The meanings of parameters in the system (2.1) are shown in Table 1. (2.1) Recovery rate of infected persons K(x) Half-saturation concentration

η(x)
Effectiveness of vaccine Remark 1 We considered that there is a vaccine coverage rate r(x) for susceptible path S, and unvaccinated susceptible persons will be injected into the infected path with a transmission rate β(x). The susceptible person who has been vaccinated enters the vaccinated compartment, suppose the effectiveness rate of the vaccine to be η(x), and if the vaccine fails, the vaccinator will also be injected into the infected path since the inoculated host has some resistance to the virus after being vaccinated. Thus this propagation process is assumed to obey a half-saturation rate (1- In addition, because R(t) does not appear in the first three equations of (2.1), we denote system (2.1) as with the initial value and boundary conditions It is sufficient to determine the dynamics of (2.1). Here, is a smooth bounded region in R n . Define a Banach space X := C( , R 3 ) with the supremum norm · and X + = C( , R 3 + ). Next, we mainly analyze the dynamic behavior of system (2.2). Let ). For any ϕ ∈ C( , R), T i is given by the following formula: and where i (i = 1, 2, 3) is the Green function associated with the operator ∂n ∂t = n in subject to the boundary condition. With [14,Section 7], T = (T 1 , T 2 , T 3 ) are compact and strongly positive. Set Then we can rewrite (2.2) as the following integral equation: For a positive and continuous function ζ (x) on , define Thus, for the local solution of (2.2), we have the following.
The proof is shown in Appendix A.
In the remainder of this section we will prove the global existence and boundedness of the solution. Consider the following equation: where D > 0 and (·), μ(·) are positive and continuous functions on . Thus we have the following.
The following theorem proves the boundedness of the model.
where (t) is the semiflow associated with the solution. Moreover, (t) is ultimately bounded.
The proof is shown in Appendix B. From what has been discussed above, we can get the following results.

Lemma 2.3
The semiflow (t) : X + → X + admits a compact and global attractor.
Proof With Theorem 2.1 we ensure the ultimate boundedness of system (2.2). Notice that the equation of (2.2) has the diffusion term, which ensures that (t) is compact. With the direct consequence in [15,Theorem 2.4.6], we can complete the proof.
In Sect. 3, we analyze the threshold dynamics of system (2.2). The global stability of disease-free equilibrium (DFE) and endemic equilibrium (EE) is analyzed by establishing the relationship between the basic reproduction number 0 and the principal eigenvalue.

Threshold dynamics 3.1 Basic reproduction number
In this section, applying the methods in [16,Sect. 3], we give the basic reproduction number 0 of (2.2). It is easy to see that (2.2) admits a disease-free equilibrium (DFE) Linearizing (2.2) at P 0 , we can get the following system: In order to discuss the basic reproduction number 0 , we will focus on the linearized equation of infected person I.
Substituting I(·, t) = e λt δ(·), we consider the following subsystem: Thus, define then the next generation operator is defined as L = -FB -1 . Set T trans as the C 0 -semigroup associated with B, then We can use the variational formula to give the solution of the eigenvalue problem as and (L) is the spectral radius of L.
Define the principal eigenvalue of (3.3) as λ 0 . Thus, we have the following result. The proof is shown in Appendix C.

Extinction of disease
In this subsection we give the proof of 0 < 1, the stability of DFE.
By the comparison principle, we can find a positive constant ϑ 1 which satisfies where ξ υ 1 is a strong positive eigenfunction associated with λ υ 1 0 . Since λ υ 1 0 < 0, we directly have This completes the proof.

Disease persistence
In this section, we prove the global asymptotic stability of the endemic equilibrium in the case of 0 > 1. First, using [10, Lemma 4.1], we get the following conclusion, which ensures that (2.2) has a positive epidemic equilibrium. The proof is shown in Appendix D.
Next, we conclude this section by proving the global stability of the endemic equilibrium.

Theorem 3.2 For 0 > 1, (2.2) admits at least one positive steady state, and we can find positive ε for any
Proof Define two sets as With Lemma 3.3, for φ 2 ∈ H 0 , we can find that x ∈ , ∀t ≥ 0, which implies I(·, t, I 0 ) > 0 and (t)H 0 ⊆ H 0 . Set here ω(φ) is an omega limit set.

Optimal control
In the previous chapters, we have focused on the disease-free equilibrium and the infectious equilibrium of infectious diseases. But if there is a sudden outbreak, we need to control the impact of the disease at a lower level as far as possible, that is, to control the number of infected people. In addition to calling for vaccinations, governments often spend more money on treatment. The mathematical language to describe this method is the optimal control problem. The main aim of this section is to develop effective strategies for controlling the spread of infectious diseases. We hope that the number of infected people does not exceed the number of susceptible and effective vaccinators.
In this section, we introduce the control strategy to (2.2) and analyze its properties. For convenience, we rewrite (x) as , the same for other parameters. To complete our research, we analyze the control variables of the model (2.2). Therefore, the control variables are given as follows.
With the development of medical technology, infected patients can be treated better. Therefore, define u = u(x, t) represents the medical intervention for infected patients. Considering that medical resources are limited, we use cuI 1+ωI for specific. Here, c is the cure rate and ω denotes the saturation constant.
From this, we give the control system of (2.2) as follows: with the boundary condition here u is measurable, other parameters are the same as in (2.2). Define an objective function t). Assume that the control set U ( × [0, T]) is convex, A 1 , A 2 are weight of each item. This objective function describes our goal to control the problem: to reduce the number of susceptible and infected people with minimal intervention costs. The value function is defined as J 0, φ(·, 0); u(·, t) .
Define a Hilbert space H = L 2 ( 2 ) and S 0 c , I 0 c , V 0 c > 0 as the initial value of (4.1), which satisfies (IOC: Initial value of Optimal Control) With the method in [19,20], we have the following lemma. Lemma 4.1 ensures the existence and boundedness of the global solution of system (4.1). Next, according to [21], we can analyze the existence of optimal control of system (4.1).

Theorem 4.1 Let the initial value be defined in (IOC).
Then there exists an optimal solution P = (S , I , V ) of the control system (4.1) corresponding to optimal control u .
Proof From the boundedness we have proved in Sect. 3, we can infer that p = inf J (u(x, t)) is finite. Thus, for P n = (S n , I n , V n ) and u n ∈ U ,we can find a sequence (P n , u n ) that is the solution to the following subsystem: with the initial condition S n (0, t) = S 0 (x), I n (0, t) = I 0 (x), V n (0, t) = V 0 (x) and the boundary condition With Lemma 4.1, for system (4.2), we can infer that ∂S n ∂t L 2 (Q) + S n L 2 (0,T;H 2 ( )) + S n (t) H 1 ( ) + S n L ∞ (Q) ≤ C, ∂I n ∂t L 2 (Q) + I n L 2 (0,T;H 2 ( )) + I n (t) H 1 ( ) + I n L ∞ (Q) ≤ C, Since H 1 ( ) is compactly imbedded in L 2 ( ), we can also get the compactness of S n , I n , V n and ∂S n ∂t , ∂I n ∂t , ∂V n ∂t . Here, by using the Arzela-Ascoli theorem [22], for the compactness we proved in Sect. 2, uniformly in L 2 ( ) with respect to a subsequence denoted by P n . In addition, with the weak convergence of S n , I n , V n (with the boundedness in system (4.2)), we have Next, we focus on the second equation of (4.2). By direct calculation we have (1r)βS n I n + (1η) αV n I n K + I n -(1r)βS I + (1η) αV I K + I = (1r)β S n I n -S I + (1η) αV n I n (K + I n )(K + I ) K + I I n S n -S + KS I n -I .
Similarly, we can discuss the first and third equation of (4.3). For the subsequence u n , u n → u weakly in L 2 (Q). With the convexity and closeness of U , hence u ∈ U , we have cu n I n 1 + ωI n → cu I 1 + ωI .
With the analysis above, we can give the conclusion that for n → ∞, P = (S , I , V ) is the optimal solution associated with the optimal function u of system (4.1), which completes the proof.
The Hamiltonian function of the control system is given as follows: H (x, t, S, I, V , u, p) (4.4) Next we give the adjoint equation for the control system (4.1) (4.5) Using the method in [23], give the partial derivative of the Hamiltonian function to u, substitute it into the optimal control solution P , Let ∂H ∂u = 0, the optimal control pair u satisfying the minimum value of the objective function min u(x,t)∈U J(u) can be expressed as u = min max p 2 (x, t)cI A 2 (1 + ωI ) , 0 , 1 .

Numerical simulation
In this section, we use numerical simulation to verify the stability of the system and the impact of controls on the development of the disease. The values of each parameter are shown in Table 2.

Stability of equilibrium
In this section, we discuss the stability of the solution of system (2.2). First of all, use the method in [28] to make difference and solve system (2.2). Then, for 0 < 1 and 0 > 1, take the data in Data 1 and Data 2, respectively. We get the simulation results obtained as follows.
Then, using Data 2, we get the numerical simulation when 0 > 1. With Theorem 3.2, we can get the uniform persistence of disease. Actually, as we can see in Fig. 2(b), with t → ∞, I(x, t) approaches a constant positive value. This is consistent with our proof in Theorem 3.2.

Influence of control parameters on disease progression
In this section, we discuss the impact of control on a disease. The main idea of this section is that we solve the optimal control problem by applying the iterative method. Then the optimal system is obtained by using the state equation and adjoint equation given in Sect. 4.
And by solving the optimal system, the optimal control strategy is obtained. Furthermore, the method in [28] is used to make difference and solve the control system and the adjoint equation. In order to control the susceptible population, the infected population, and the vaccinated population, targeted treatment of patients is taken as control, and the impact of targeted treatment on the susceptible population, the infected population, and the density of the vaccinated population in a long-term state is considered. Finally, through numerical simulation, the actual situation of each path in the original system (2.2) and the control system (4.1) is compared.
First, we define the objective function as follows. Let the objective function corresponding to the control system (4.1) be as follows: where A 1 = 0.4, A 2 = 0.5 [29]. The values of other parameters are the same as in Sect. 5.1.
The numerical simulation results are as follows. In Fig. 3, for 0 < 1, under controlled conditions, the duration of the disease is shorter and the time of extinction is earlier than without control. At the same time, the overall density of susceptible and vaccinated people before reaching the stabilization point was also higher in the controlled condition than without control. In addition, according to Fig. 4, the control intensity reached the maximum in the early stage and gradually decreased with the weakening of the disease scale, reaching zero value when the disease disappeared. The optimal control when 0 < 1 Figure 5 Infectious path S, I, V without and with control ( 0 > 1) Figure 6 The optimal control when 0 > 1 As shown in Fig. 5, for 0 > 1, when control exists, the scale of the disease reaches a minimum earlier than without control, and when the disease eventually becomes endemic, the total scale of the disease is lower than without control. In addition, when the disease reaches a stable state, the density of susceptible and vaccinated persons is higher in the controlled condition than without control. Furthermore, by Fig. 6, when the disease is in its initial state of development, control rises, and when the disease reaches equilibrium and becomes endemic, control is maintained at a stable value along with the duration of the disease.

Conclusion and discussion
In this paper, a kind of SIVR infectious disease model including vaccine immunity and vaccine effectiveness is considered. The optimal control theory is applied to the study of the model, and the threshold dynamics and optimal control of the model are discussed. Firstly, we prove the well-posedness of the model, which provides a theoretical basis for the following discussion. Secondly, we give the basic reproduction number 0 to analyze the dynamic behavior of the disease threshold. In addition, the Hamiltonian function and adjoint equation of the optimal control problem is given. Finally, the stability of the system solution is verified by numerical simulation and the number of infections can be reduced as much as possible, while the cost is reduced under the treatment control. In this paper, the parameters are assumed to be accurate; in fact, due to various uncertainties, each parameter may be inaccurate or random. At the same time, according to the changes in the parameters, it can be seen that the vaccination rate and the effective rate of the vaccine also have a certain impact on the control (see Fig. 7(a), (b)). In addition, because the nearoptimal control is more flexible, it can adapt to different degrees of model uncertainty. Therefore, in future work, the near-optimal control problem of the epidemic model can be further studied by considering the influence of random parameters, noise, the vaccination rate, and the efficiency rate of the vaccine as the control parameter.